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Summary. The theory of cointegration has been the leading theory in econometrics with powerful 
applications to macroeconomics during the last decades. On the other hand the theory of phase syn- 
chronization for weakly coupled complex oscillators has been one of the leading theories in physics 
for some 15 years with many applications to different areas of science. For example in neuroscience 
phase synchronization is regarded as essential for functional coupling of different brain regions. In 
an abstract sense both theories describe the dynamic fluctuation around some equilibrium. In this 
paper we point out that there exists a very close connection between both theories. We show that 
there exists a system of stochastic difference equations which can be viewed as close to the Ku- 
ramoto equations, whose solution is a cointegrated system. As one consequence the rich theory on 
statistical inference for cointegrated systems can immediately be applied for statistical inference on 
phase synchronization based on empirical data. This includes tests for phase synchronization, tests 
for unidirectional coupling and the identification of the equilibrium from data including phase shifts. 
We give an example where an unidirectionally coupled Rossler-Lorenz system is identified with the 
methods from cointegration. The methods from cointegration may also be used to investigate phase 
synchronization in complex networks. On the contrary there are many interesting results on phase 
synchronization which may inspire new research on cointegration. 

Keywords. Cointegration, phase synchronization, weakly coupled oscillators; driver response rela- 
tionship; Rossler-Lorenz system. 



Address for correspondence: Institute of Applied Mathematics, University of Heidelberg, Im Neuenheimer 
Feld 294, D-69120 Heidelberg, Germany (e-mail: dahlhaus@statlab.uni-heidelberg.de). 



Contents 

1 Introduction 

2 Phase synchronization and the Kuramoto model for interacting oscillators 

2.1 Phase synchronization of weakly coupled oscillators 

2.2 Cointegration as a stochastic Kuramoto-type model 

3 The concept of cointegration 

3.1 Cointegration and the Johansen-Granger representation theorem 

3.2 Cointegration Testing and Model Fitting 

4 The cointegration approach to phase synchronization 

5 Example: A coupled Rossler-Lorenz system 

6 Conclusion and outlook 

7 Appendix 1 : VEC - state oscillators 

8 Appendix 2: Computational and modeling aspects 

8.1 Computational aspects 

8.2 Modeling aspects 







1 Introduction 



Phase Synchronization has a long history dating back to 1 665 where the mathematician and 
physicist C. Huygens discovered some antiphase synchronization of two pendulum clocks 
suspended close to each other on the same wooden beam. From that time on, the phe- 
nomenon got more and more into the focus of scientists. During the last decades the be- 
havior of two or several interacting oscillators has been intensively studied in the physics 
literature - in particular in the context of nonlinear dynamical systems (cf. [59, l42l 1701 19118): 
- see also the monographs [63l [53]. In contrast to other types of synchronization, phase 
synchronization purely depends on the phases of self-sustained weakly-coupled oscillators 
while the amplitudes may even be empirically uncorrelated. Weak coupling means that the 
phases of all oscillators may be subject to individual disturbances and the whole system 
adjusts itself afterwards. Thus phase synchronization is rather regarded as a complex dy- 
namical process than a fixed state. 

The phenomena of phase synchronization has been found experimentally in many fields, 
e.g. electrical circuits (66] [6], lasers [14], electrochemistry [401 [4~T|, biological systems [T8l 
75], population dynamics [7], between El Nino and the Indian monsoon [50] and even tennis 
[57). Since phase synchronization purely depends on the phases of the oscillators it can 
also be detected between oscillators which are of a different type. A practical example of 
this is the phase synchronization of the cardiac and respiratory system 071 [73] or between 
maternal breathing and the fetal-maternal heart rate coordination [76). 

In neuroscience phase synchronization is regarded as essential for functional coupling of 
different brain regions. Single neuronal activity, phase coupled to an ensemble of oscillatory 
neuronal activity, enables the cells to transmit their information content long range across 
different cortical areas. To cite from the abstract of Fell and Axmacher [22]: "In recent years, 
studies ranging from single-unit recordings in animals to EEG and MEG studies in humans 
have demonstrated the pivotal role of phase synchronization in memory processes. Phase 
synchronization - here referring to the synchronization of oscillatory phases between differ- 
ent brain regions - supports both working memory and long-term memory and acts by facil- 
itating neural communication and by promoting neural plasticity." [10] discuss mechanisms 
of gamma oscillations in the hippocampus and the functional role of the synchronization of 
such oscillations in several key hippocampal operations, including cell grouping, dynamic 
routing, and memory. In [79] synchronization in the gamma band is investigated, and it is 
discovered that the mutual influence among neuronal groups depends on the phase relation 
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between rhythmic activities within the groups. Furthermore, the pattern of synchronization 
is related to the pattern of neuronal interactions. Further reviews on phase synchronization 
in neuroscience are (1911771 [13) . 

To study the synchronization of a larger population of interacting units theoretically - for 
example fireflies flashing at the same time or crickets chirping at the same time - Winfree 
[78] studied the nonlinear dynamics of a large population of weakly coupled oscillators with 
intrinsic frequencies that were distributed according to some prescribed probability distribu- 
tion. He also ignored the amplitude and considered phase oscillators and worked with a 
mean field model. Kuramoto [44, 45] introduced a more sound mathematical model to de- 
scribe this phenomenology leading to a theoretical treatment of the mean field approach. He 
also studied the limit behavior when the number of oscillators tends to infinity. This model 
has been used since then in different forms to discuss theoretically phase synchronization 
in a population of weakly coupled oscillators. For example the onset of synchronization is 
discussed in this framework [74]. 

Although some statistical methods are used to analyze phase synchronization (e.g. the 
spectral coherence) there hardly exist any stochastic models for synchronized phase pro- 
cesses. Such stochastic phase models would be an important task to estimate the dynam- 
ics of the phase processes from observed data. Specific outcomes could be tests for phase 
synchronization, tests for unidirectional coupling and the identification of the equilibrium from 
data including phase shifts. 

In this paper we argue that the famous theory of cointegration provides also a good 
stochastic model for phase synchronization and therefore is a good framework for inves- 
tigating these questions. The trick is to use cointegration not directly as a model for the 
oscillators but as a model for the phase processes driving the oscillators. Using the phase 
processes as the key element for describing oscillators has been standard in physics for 
many years. On the other hand this approach has never been used to our knowledge in 
statistics or econometrics to model oscillators. Here oscillators have always been investi- 
gated by means of Fourier transforms (spectral methods) or wavelet transforms. Further- 
more, statisticians have always tried to build models directly for the oscillating processes 
instead for the phase processes. 

The theory of cointegration may be used in phase synchronization both for theoretical 
studies and for identifying unknown systems based on empirical data. The link may also 
be of importance in the other direction in that the known results on phase synchronization 
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may lead to new insights or stimulate research in cointegration. The common ground of both 
theories is that two or several processes fluctuate randomly around some equilibrium. 

The concept of cointegration was introduced by Granger 1981 [24] for the joint modeling 
of several macroeconomic variables over time. In 2003 he received the Nobel Prize for his 
discovery jointly with Robert Engle. The aim of the theory is to model common stochastic 
trends, for example of income and consumption where the short-run dynamics is affected 
by random disturbances and the long-run dynamics is restricted by economic equilibrium 
relationships. Other examples are exchange rates and price levels, short and long-term 
interest rates, and the prices on spot and future markets. 

Since its introduction the theory of cointegration has been developed by many researchers 
and it has become the leading theory in econometrics with many applications in macroeco- 
nomics and beyond - for references see the monographs (5j |35] HH @9] 36] among many 
others. |61]|39] have developed the theory for continous time diffusion processes. 

We now give a heuristic and elementary argument why the setting of cointegration is of 
importance for phase synchronization: Suppose we observe yf' = AiCos(<ffl) + ef 3 for 
i = 1, 2 and t = 1,...,T with phase processes <ffi (alternatively we may just have the 
phases 4>t from some oscillators calculated by means of the Hilbert transform or some 
other method). A naive model for a random phase could be <j> t = wt + 0o + S t with some 
random error 6 t which could be a white noise process or some stationary process. This 
model had the drawback that it fluctuates around the deterministic linear phase ojt + <fi and 
would therefore only be adequate if some external force or constraint would keep the phase 
close to this linear phase. Instead a better model for most situations is the corresponding 
model for the phase increments A<f) t = <p t - <^_i = u + 5 t where 5 t is a stationary process 
with mean zero and positive correlation. In the simplest case where 6 t is iid Gaussian this 
means that <j> t = ut + + J2l=i &s has the same distribution as a Brownian motion with 
drift. In this case var(0() ~ t which means that the phase may depart substantially from 
4>t = cot + 4>o for large t. This feature remains if the process 8 t is stationary with mean 0. 
Such a process 4>t is called an integrated process. 

If we look at two processes with synchronized phases then both phases may evolve like a 
Brownian motion with drift - however the difference 4>y - (j)^ should stay relatively small; in 
particular it should not explode. A proper model therefore is that these differences follow a 
stationary process. This is exactly the concept of cointegration: both phase-processes are 
integrated but the difference remains stationary. The heuristics presented here is formalized 
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in Appendix 1 . 

The paper is organized as follows. Section [2] provides a brief introduction to phase syn- 
chronization with focus on testing for phase synchronization and the Kuramoto model for 
interacting oscillators. We point out that a cointegrated system arises as the solution of a 
system of stochastic difference equations which can be viewed as close to the Kuramoto 
equations. In Section [3] we give a brief introduction to the theory of cointegration and review 
some tests for cointegration. Section @] describes the use of cointegration for the statis- 
tical analysis of phase synchronization. In particular we give our definition of stochastic 
phase synchronization. In Section [5] we present a simulation showing that the proposed 
method works for a coupled Rossler-Lorenz system. Conclusions and an outlook are given 
in Section [6]followed by the definition of VEC-state oscillators in Appendix 1 and some com- 
putational and modeling aspects in Appendix 2. 

The article wants to address physicists, econometricians and statisticians, as well as 
scientists who just want to apply the method. 

A quick reading guide for physicists may be to start with the discussion of the Kuramoto- 
model in Section I2T21 to see how ®, (|9j, (T4}, and (T5) are related to the Kuramoto-model 
and to proceed with Sections Hand [5] 

A quick reading guide for econometricians and statisticians familiar with cointegration 
may be to start with the model for a VEC-state oscillator in Appendix 1 , continue with the 
example in Section |5]and to pass on to the Kuramoto-model in Section |2]2]and the definition 
and discussion in Section HI 

Scientists mainly interested in the method may start directly with the example in Sec- 
tion |5]to get a first impression on cointegration and on the impact of cointegration for phase 
synchronization. 

2 Phase synchronization and the Kuramoto model for interact- 
ing oscillators 

2.1 Phase synchronization of weakly coupled oscillators 

The starting point of phase synchronization is the definition of the phase 4>t of an oscilla- 
tor. Traditionally phase synchronization has been mainly considered for chaotic oscillators 
where the phase often is defined via the Hilbert transform of the signal. If the projection 
of the attractor on some plane has only one rotation center then the phase can be defined 
alternatively by using a Poincare map or a certain rotation angle (64). Alternatively, one may 
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look at stochastic systems, e.g. of the form 

y t = a t cos(0 t ) + e t , t£Z (1) 

or even at systems like 

Vt = a t f((j)t) + bt + e u teZ (2) 

where / is a 2ir periodic real valued function representing the oscillation pattern, and a t and 
b t are the amplitude and the baseline respectively [T2]. The latter model may for example 
be used for ECG signals. 

Based on discrete data of the oscillator, the phases <p t have mainly been determined via 
the Hilbert transform [70]. The Hilbert phase of a signal y t is defined through the analytic 
signal Q given by Ct = Vt + W? = a t exp(i4> t ). The imaginary part of ( t is obtained using the 
Hilbert transform which is defined through 



7T J t — S 

with P being the Cauchy principal value. a t is the amplitude and the phase 4> t is computed 
through 

(j) t = arctan— . 
yt 

Other methods for phase estimation are based on the Wavelet transform (26), on tne P eri - 
odogram [29j, on a local periodogram [58] or on a general state-space model with particle 
filtering [T2]. The last method has in particular been used in combination with the model ©. 

In the next step two weakly coupled oscillators with phases <j>jp and 4>f' are called m -. n 
phase synchronized if there exists some constant <I> (phase shift) and some small 5 > such 
that 

| (m4 l) ~ n4 2) ~ mod 2vr | < 5 (3) 

holds for all t. m and n are integers which are usually known in a practical application. The 
idea behind this definition is that the (rescaled) phases "stay together" and do not move 
arbitrarily far from each other. A standard measure for testing for phase synchronization is 
the phase synchronization index (51] [67] defined through 
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-f exp { 1 i 171 ^ - n( rf 2) ) } 

The synchronization index is an estimate of the corresponding theoretical quantity called 
mean phase coherence. It is normalized between and 1 where values close to one in- 
dicate phase synchronization. A value close to 1 is obtained for an almost constant phase 



difference of —n<ffi | which has also been observed between two nonidentical chaotic 
oscillators (70]. R 2 has often been used for testing phase synchronization and there exist 
several articles where the distribution of R 2 or other test statistics is approximated by dif- 
ferent methods - e.g. with surrogate data. The most advanced approach is (72) where the 
distribution of R 2 under the null hypothesis R 2 = is approximated by using two indepen- 
dent stochastic processes. In [2] different tests are presented by using a two sample test 
setup and by utilizing insights and methods from directional statistics and bootstrap theory. 

Comprehensive reviews and comparisons of various phase synchronization measures 
and tests can be found, in (67), 0. an d |72) . 

Another important aspect of phase synchronization is directional dependency and there 
exist approaches based on information theory with surrogate data [55] and mutual pre- 
dictability [69] to detect the direction of coupling. A practical example of a phase syn- 
chronized system with unidirectional coupling is the cardiac and respiratory system. (69) 
provided evidence that the respiration is driving the cardiovascular system. We mention that 
many other directional measures exist which are associated with other types of synchroniza- 
tion (cf. [56] and the references given there). 

Despite of these results a good stochastic model for the phase processes themselves in 
which questions like synchronization and unidirectional coupling can be investigated is still 
missing. We show in this paper that cointegrated systems provide such a framework. 

2.2 Cointegration as a stochastic Kuramoto-type model 

As mentioned in the introduction Kuramoto HU |45) introduced a mathematical model to 
describe the synchronization of a larger population of interacting units. A simple version of 
the Kuramoto model for the phases of d oscillators is 

^ = w ' + f E sin ~ (i = h- ■-,<!). (4) 

3=1 

Allefeld and Kurths [1] discuss the following stochastic generalization of the model 

d 

^=u i + J2k ij wx(^-<f,^)+TH (i = l,...,d) (5) 

3=1 

where the rji are taken to be independent Gaussian white noise variables. 

There is a striking similarity of this model to cointegration: For <f>W — </>W small we can 
approximate sin(^) - <j>^>) by ^C?) - leading to the discretized approximation for the 
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vector fa = [4>f\...,ct>f ) )',u= (wi,...,^)' 

= w + n&_i +7/t (6) 

with n = if-diag(X)^=ifcii,- • • ,Ylj=i k dj) where K = , d . This is exactly the 

representation © from below with p = 1 and w = n, i.e. the cointegrated system fTO), © 
with p = 1 is exactly the solution of the stochastic difference equations ©. We mention that 
cointegration requires the matrix n to be of reduced rank which is fulfilled in this case since 
the columns sum up to 0. 

Being a bit more specific at this point gives us a deeper insight into the relevance of 
cointegration theory for phase synchronization in networks: The above n can be written in 
the form n = aft with d x r - matrices a and p of full rank r where (in this case) r = d — 1, 
namely 

(I ••• -1\ 
1 ••• 0-1 

\0 ■■■ 1 -1/ 

and a consisting of the first d - 1 columns of n (only for this specific /3). The ftfa = 
are obviously the d - 1 phase synchronization relations ffl = ■■■ = The matrix a is 
sometimes called adjustment or loading matrix. 

In the following we apply the concept of cointegration to phase synchronization. For the 
model © this means testing for cointegration (phase synchronization) via determining (test- 
ing) r = rankll. This results in the reduced rank representation U = a/3' with dxr - matrices 
a and /3 of full rank r. After the parameters a and (3 have been estimated ft fa are the r 
cointegration (phase synchronization) relations and a gives the intensities with which devia- 
tions from the equilibrium lead to corrections. In particular a can be tested for unidirectional 
coupling. 

In addition it is very interesting to see that the vector u splits under cointegration to a 
linear drift term oj (basic frequency) and a constant term in the equilibrium (3 (mean phase 
shift) transforming © finally to 

Afa = ujQ + a (ftfa-i - Pa) + Vt 

meaning that the "true phase synchronization relations" are ftfa-i - p = 0. Of course the 
corresponding Kuramoto model © and © should then contain sin(0^ -Pq) instead 
Of sin(c/.(^ - 0^). 
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In the example of Section [5] we apply this analysis with the more general cointegration 
model Ci3 in the case d = 2 to data coming from a coupled Rossler-Lorenz system. 

There is some connection to the work [4] where also the linearized dynamics of the Ku- 
ramoto model is investigated (however in a deterministic framework). As in that work the 
phase synchronization relations come out as the eigenvectors of some eigenvalue problem 
with the strength of the synchronization being determined by the magnitude of the eigenval- 
ues. Due to the stochastic structure the eigenvalue problem however is different. 

Summarizing: in the case of phase synchronization, where (say) <j>^ - <t& is small, the 
model © is almost equivalent to a stochastic Kuramoto model; in the case of non phase 
synchronization the model © (or better its generalization ©) is also a proper stochastic 
model for the phase processes (in this case rankn = or d). As a consequence the model 
© can be fitted in both situations (phase synchronization / non phase synchronization ) 
to the data and the hypothesis of phase synchronization can be tested on the matrix n - 
e.g. by looking at the fitted model. In addition one is able to identify the parameters of 
the "Kuramoto-type" model from real data, and to conclude to the phase synchronization 
relations, undirectional coupling etc. 

Under the hypothesis of phase synchronization the difference between the two models 
©, © (resulting from the replacement of sin(<^— <f>®) with ^CO— <^W) is small but of course 
not negligible. In the Kuramoto model the phase differences are employed modulo 2ir while 
in the cointegration model the ordinary differences between the phases are employed. In 
practise this means for example that two processes whose phase difference is close to zero 
at the beginning but moves at a certain time point (quickly) to 2ir may be phase synchronized 
in the classical sense but not phase synchronized in the cointegration based definition given 
below (this happens exactly in the example of Section [5]for e = 9.6, 10.2 - see FigureE}. 

3 The concept of cointegration 

3.1 Cointegration and the Johansen-Granger representation theorem 

We now briefly review a small part of the theory of cointegration (for more details see e.g 
the monographs [35l|49][36] - the last one containing a more applied view). As it is usual in 
many papers we restrict ourselves to vector autoregressive (VAR-) processes 

X t =A 1 X t -i + ... + A k X t - v + n + TH, (7) 
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where the coefficients Ai are d x d-matrices, \i = (m, . . . , /j,d)' and the errors are inde- 
pendent identically distributed (iid) random variables with mean and positive definite co- 
variance matrix 9, n . Let A(z) •— h - A\z A p z p be the characteristic polynomial. We 

restrict ourselves to VAR-processes with roots outside the unit circle or equal to 1, i.e. where 
det (A(z)) = implies \z\ > 1 or z = 1. 

We call a process integrated (of order 1) (1(1) for short) if X t is not stationary and the 
first order difference AX t = X t - X t -i is stationary. We call a d-dimensional process X t 
cointegrated if each univariate series is integrated and some linear combination j3'X t (with 
an r x d-matrix /3 of rank r and 0<r<<i) is stationary (these definitions simplify the issue - 
for thorough definitions we refer to the above monographs). /3 is known as the cointegrating 
vector. It specifies the long-term relationship between the involved univariate series. The 
most relevant example for phase synchronization of 2 oscillators is d = 2 and (3 = (1, -1)': 
(3'X t fluctuates stationary around a long run mean - i.e. it is bounded in probability. 

With some straightforward calculations we can transform the above representation to 

p-i 

AX t = n X t _i + Y, T i Ax t-i + V + Vt (8) 

i=l 

where n := -(I d -A 1 A p ) andTj := -(A i+1 -\ h^4 p )fori = l,...,p-l. If the process 

X t is not stationary (in particular if it is cointegrated) then det(n) = (-l) d det (-4(1)) = 0. 
Thus n is singular with rank r < d and it can be decomposed to n = a/3' with dxr - matrices 
a and (3 of full rank r. 

Since the process is 1(1), the first order difference AX t is stationary, and the represen- 
tation |8]) implies that also UX t -i = a/3'X t -i is stationary. Multiplying this with (a'a^a' 
implies that f3'X t -i is stationary - i.e. each component of p'X t -i is a cointegrating relation. 
The corresponding representation 

p-i 

AX t = a p'X t „i + r* AI M + n + rit . (9) 

i=i 

is called vector error correction model (VEC-model or VEC-representation). It shows that 
whenever the process moves away from the equilibrium ft'Xt-i = fio (for the definition of (3 
see below) it is pulled back to the equilibrium with the forces a (see the end of Section [5]for 
an illustrative example). 

The Johansen-Granger representation theorem [20, 33] now gives another useful repre- 
sentation of the process. Under the above conditions the process has the moving average 
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representation 

t 

X t = C'Y,Vi + Cfit + C*(L)( Vt + f i)+XZ (10) 
i=i 

where 

C = pAa'jTpX a' L (11) 

with r := I p - J2 P i=l IV C*(L)(r, t +n) = 22jLo Cj(m-j + ^) is a stationary process and X* 
contains initial values with /3'Xq = 0. 

aj_ denotes an orthogonal complement of a i.e. a ± is an d x (d - r)-matrix of rank d — r 
with a' ± a = (the same for f3±_). 

The representation (TO) decomposes the process X t into nonstationary 1(1) components 
and stationary components. It shows that the d-dimensional process is driven by d-r 1(1) 
components and r stationary components. The first term on the right hand side consists 
of d random walks Yh=i m which are multiplied by a matrix of rank d-r denoted by C. 
Thus, there are actually d-r stochastic trends driving the system. On the contrary the 
representation (9) (or even better the representation {T3) below) shows how the process is 
pulled back to the equilibrium if deviations occur. This is illustrated in more detail with the 
example at the end of Section [5) 

Of special importance for phase synchronization is a decomposition fi = Toj - a(3 where 
w is a trend-term and ^ is a constant belonging to the error correction equation: Taking 
expectations in ([TO]) we obtain because of Brj t = 

BX t = Cut + C*(L)fi + EX£ (12) 

and therefore w : = EAX t = C\i. Since AX t is stationary we obtain from (9]) 

p-i 

aE(A-i) =E(AX t ) -^r i E(AX t _ i ) - fx = [TC-I d ]fJt. 

i=l 

Let a := a(a'a) -1 . Since a' a = I r we have E(/3'X t -i) = a'[TC - I d ]n =: Po- 

Obviously we have \x = TCfi - [TC -I d ]n = Tuj - a/3 Q . Therefore we obtain the modified 
VEC-representation 

p-i 

(AX t -uj ) = a - /3o) + T i ( AX t-i ~ w o) + Vt- (1 3) 

i=l 

This means that the "true" cointegration relation describing the equilibrium is fi'X t -i-Po = 0. 
If the process deviates at time t — l from this relation it is pulled in the next step with force 
a back towards this equilibrium. (Q7jJ and {T2]) show that u = Cfj, is the drift-vector of 
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the process, that is the process X t has a deterministic trend uj t. The intercept terms are 
contained in Xq in (pfO). 

Below we use this model as a stochastic model for phase synchronization. In the case 
d = 2 the situation is even more intuitive: If the system is cointegrated (i.e. r = 1) we only 
have a one-dimensional cointegration relation p'X t -i - (3 = and the drift vector u has 
the same direction as the random walk part CYX=i m (namely f3 ± ). We discuss this simpler 
and more intuitive case in Section [5] with a specific example. 

To keep the situation simple we also restrict ourselves to the case p = 2. In this case the 
VEC-representation takes the form (now replacing X t by <j> t ) 

A0{ 1} = a 1 {4 1 2 1 -p 2 4 2 \)+ 11 A4 1 2 1 + 8 1 A ( l ) ^ 1 +^ + 4 1) , (14) 
A4 2) = a 2 (4\-^4 2 2 1 )+^A4% + S 2 A4 1 } 1 + f , 2 + 4 2 \ (15) 

where we have assumed for the cointegrating vector /3 = (1, -f3 2 )' for identifiability. In that 
case also a\ and a 2 are identifiable. 

If the model is used for phase synchronization j3 2 usually is known (in the 1 : 1 synchro- 
nization case we have (3 2 = 1; /3 then is the mean phase difference in the equilibrium). We 
expect that in many cases the reduced system with Si = 5 2 = will suffice to describe the 
joint phase dynamics (this is for example confirmed by parameter testing in the example of 
Section 5 below). However, the autoregressive part is needed since successive increments 
usually are correlated. 

It can be seen that the error correction coefficients a = (ai,a 2 )' influence the speed of 
correction. The Johansen-Granger representation theorem implies «i < 0, a 2 > or both if 
4>t and <j>y are cointegrated. Note that this allows one to determine whether there is an 
unidirectional driver-response relationship between the processes by testing if one of the 
error correction coefficients is zero. 

3.2 Cointegration Testing and Model Fitting 

In order to test for phase synchronization and for unidirectional driver-response relationships 
we shall apply tests for cointegration. Two cases need to be distinguished: First the case 
where the cointegrating vector /3 is known and we only want to test for a specific cointegrating 
relationship. This may be the standard case for phase synchronization of (say) two oscilla- 
tors where the m-.n synchronization relation as in (3) is often clear from prior knowledge or 
eye-inspection. 
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The other case is the case of unknown (3. For example in systems of higher dimension it 
may be clear that the whole system is "somehow" 1 : 1 synchronized but it usually is not clear 
at all how the synchronization "propagates" through the network. In that case the task is to 
test for the rank of n in ® and to determine the factorization n = a ft as well as the vector 
(3 leading to the phase synchronization relations, the mean phase shift and the matrix a 
representing the strength and the direction of coupling. Some examples and a discussion of 
the meaning of the rank(n) can be found in the next section. 

We start with the case of unknown f} leading to the likelihood theory of Johansen (35) 
and in particular to his likelihood ratio test. Additional information on this case is provided in 
Section m We then discuss the augmented Dickey-Fuller unit root test f1"5l[1"6"] . 

Johansen's likelihood theory and rank test 

The likelihood theory for cointegrated systems is well developed - the most famous result 
being Johansen's likelihood ratio test (LR-test) for the cointegration rank - i.e for the "deter- 
mination" of r = rank(n) in ® and the corresponding decomposition 11 = a p' with a and /3 
being dxr - matrices of rank r. Thus a major advantage of the Johansen procedure is that 
it can simultaneously identify multiple cointegration relations in multivariate time series. The 
situation is much more challenging than for the classical likelihood ratio test since the inte- 
grated processes require a different asymptotic theory leading for example to a nonstandard 
limit distribution of the likelihood ratio test. 

We now briefly sketch part of the results. Since the results can hardly be summarized in a 
few lines we strongly recommend reading for example the applied presentation in Chapters 7 
and 8 of [36]. When reading this presentation the reader should keep in mind that inference 
is based on the model ®,((9j and that n (leading to the deterministic drift term and the phase 
lag as described above) plays an important role. The deterministic term m discussed there 
(called drift term and leading to a quadratic trend) can be ignored (at least for "typical" phase 
synchronization problems). Thus the case to be considered is Case 3 in [36]. 

Suppose now we want to calculate the maximum likelihood estimates for a and (3 in the 
system ®. In the Gaussian case the maximum likelihood estimates are essentially the same 
as least squares estimates. In the first step // and the AI t _i are removed by regressing 
them on AX t and X t -\. If we denote the residuals by AX t and X f _i we then have the 
"concentrated model" 

AX t = a P'X^i + error. 
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Now, the optimal value is determined for given /3, plugged into this equation, and finally 
the optimal /3 is determined leading to the maximum of the likelihood. This procedure can 
be viewed as finding the canonical correlations between AX t and X t -i (c.f. [36], Chapter 
8). The resulting squared canonical correlations are denoted by Ai > . . . > A^ > 0. This 
means A\ is the squared canonical correlation of that linear combination /3[X t -i (equilibrium 
relation) having the highest correlation with the increments AX t and therefore the highest 
error correcting influence on the process (phases). Thus Xj belongs to the cointegration 
(phase synchronization) relation /3J-X f _i. If a A^ is zero (close to zero based on data) it does 
not belong to a cointegration (phase synchronization) relation. Therefore the task is to split 
the Ai > . . . > Xd into those Ai, . . . , X r which are different from zero and which belong to 
a cointegration (phase synchronization) relation and those A r +i, . . . , X d which are zero (i.e. 
where the corresponding unknown theoretical values are zero). With view on (fTTJ) and the 
discussion below, d - r is the number of driving "forces" of the trend. In particular r is the 
rank of the matrix n we are looking for. 

A formal way to determine r = rank(n) based on the values of Ai, . . . , A^ is the likelihood 
ratio test - e.g. the so-called trace-test where the null hypothesis H(r) = {rank(n) < r} is 
tested in H(d) -.= {rank(n) < d} recursively for r = 0, . . . , d - 1 (this recursion is called "top 
->■ bottom" procedure). The test statistic is 

d 

r d _ r := LR(H(r)\H(d)) = -T Ml - A*) 

i=r+l 

(called LR-test below) and the hypothesis is rejected at a significance level a if r d _ r > 
C a (d - r) where C a (d - r) is from Appendix A of [36] (Case 3; p—r = d—r). The procedure 
stops if the test accepts the hypothesis for some r. 

Special attention is required for the case r = d. This means that there are no stochastic 
trends and the system in ® is stationary. In particular this case does not mean that there 
are r cointegration relations - see the discussion in the next section. 

Augmented Dickey-Fuller Unit Root Test 

The Augmented Dickey-Fuller (ADF) test allows one to test for the presence of a unit root in 

a univariate time series Y t . It is based on the regression 

p 

AY t = a + bt + cyt-i + cjAY^j + e t . 

The constant a and the time trend bt are only included if required. For a ^ and 6 = 
the ADF-test tests the null hypothesis c = integrated process with a deterministic trend 
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against the alternative c < stationarity with a mean. The test statistic is the usual t-ratio 
given by 

SE(c) 

with c being the ordinary least squares estimator of c and SE(c) being its estimated standard 
error. The corresponding quantiles can be found in [28] where in our case the standard 
normal distribution must be used. In case the potential cointegration relation f3'X t is known 
the ADF-test provides a simple method for cointegration testing. We then apply the ADF-test 
to the residuals Y t = j3'X t -\. Due to a / it is not necessary to include fa. 

A disadvantage of the ADF-test is the implicit common factor restriction which is imposed 
when the ADF-test is used. The test looses power if the restriction is not satisfied |43] . 

Other Tests 

An alternative to the ADF-test is the Wald test which tests for cointegration via the sig- 
nificance of the adjustment coefficients a (32). Again, special critical values are required 
[32] since the asymptotic distribution of the Wald statistic is also non-standard. In (17) a 
device has been constructed which guarantees that the Wald test has an asymptotic x 2 - 
distribution under a general set of conditions. It has been shown that the Wald test is supe- 
rior/inferior to the ADF-test if the innovations rjf\ rfp in the VEC-model are strongly corre- 
lated/uncorrelated respectively [32], p. 1005). In our experience the low dimensional model 
{14}, (T5]) is already sufficient for modeling typical phase processes and it typically leads to 
almost uncorrelated innovations (see Section ©. For this reason we prefer the ADF-test in 
comparison to the Wald test. 

Several other cointegration tests have been proposed in the literature. A method which 
may also be useful in our framework is based on the single-equation conditional error cor- 
rection model. The basic idea is to rewrite the VEC-model (T4},|[T5]) as a conditional model 
given 4>f^) and a marginal model (for Under the assumption of weak exogene- 
ity a cointegration test which is solely based on the conditional model can be constructed 
|34l|30l|80]. The problem of this approach is the fact that weak exogeneity is equivalent to 
a 2 = 0. Therefore, the weak exogeneity assumption implies that only unidirectional error 
correction relationships can be considered which is a strong restriction. 

The ADF-test described in Section [3~2l can also be applied when f3 is unknown. Then, 
/3 is estimated through a regression followed by the computation of the ADF test statistic. 
Due to the estimated /3 other critical values need to be used. This approach is known as the 
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Engle-Granger two-step procedure [20]. Another important test is the Phillips-Ouliaris test 
based on residuals [62). In [71] Lagrange multiplier tests are constructed for the situation 
where a time trend is present. 

We finally mention that the ADF-test can also be used to test the hypothesis integrated 
process vs. trend stationarity. An alternative is to use the KPSS-test [46] where the hypoth- 
esis are exchanged (see the discussion in the next section and in Appendix 2). 

4 The cointegration approach to phase synchronization 

The cointegration model for phase processes 

We now use the concept of cointegration for a stochastic definition of phase synchronization. 
The key idea is that in a stochastic context phase synchronization can be described in terms 
of stationarity of the phase-differences, i.e. by cointegration. That is we propose to replace 
the fixed deterministic bound in © by a stochastic bound: The difference may even get large 
(with small probability) but the fact that the difference is stationary will always force it back 
to the equilibrium. 

Before we give the definition we stress that here we only want to discuss models for sig- 
nals whose phase increments can be regarded as stationary - that is where the unwrapped 
phases are integrated processes with a deterministic linear trend (the stochastic part of the 
integrated process is sometimes called "stochastic trend"). The case where the determinis- 
tic trend is zero will hardly occur in practice; the case where the stochastic trend is zero (i.e. 
the process is trend stationary) may occur in specific cases - see the discussion below. By 
integrated we always mean l(1)-integrated in terms of the cointegration literature. 

Definition 1 (Stochastic Phase Synchronization), d oscillators with phase processes 
4> t = {4>t\- ■ ■ > are called stochastically phase synchronized of order r with 1 < r < 
d — 1, if all processes are integrated and there exist r linearly independent vectors such 
that the & <fi t are stationary (i. e. <p t is cointegrated with rank r). In the case of the VEC-model 
® this means that rank(n) = r andU = a/3' as in {g|) with = . . . , f3 r ). The p'M are up 
to a constant and up to some non-identifiability the phase synchronization relations. 

One may also use the term stochastically ji-phase synchronized of order r if the pro- 
cesses are cointegrated of rank r with cointegrating matrix/vector /3. 

We now briefly discuss the different orders r of phase synchronization which correspond 
to the rank of the matrix n in ®: 
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r = : If rank(n) = there do not exist any cointegration (equilibrium) relations and the 
phases are integrated processes with a deterministic drift term (which for phase processes 
usually is different from 0) and a stochastic drift term (meaning that the processes are inte- 
grated). Since the vector of phase increments follows a VAR-process as in ® with n = 
there still is some stochastic dependence of the stochastic drifts terms but this does not lead 
to synchronization of the phases. 

0<r<d : In this case we are having r phase synchronization relations, meaning that r is a 
measure for the degree of synchronization. We give some examples: 

- if r = 1 we may have e.g. the relation ffi = (f>f> or for example the relation 
Y^j=i ft4>i ~ Po = (this form may hardly occur in practise); 

- if r = 2 we may have two (linearly) independent relations such as 4>t = 4>t anc ' 
(j^f 1 = <pp, or two generalized relations as above; 

- if r = d — 1 we have for example the relations 4>t = ■ ■ ■ = 4>f' '■ 

Recall that this means that the d-dimensional process is stochastically fluctuating around 
these r phase synchronization relations with some error correction mechanism as given in 
(fT3) . The phase synchronization relations are not uniquely determined. This is reflected by 
the equation n = aj3' = o^' -1 ^)' with any regular r x r -matrix £ where now the columns 
of are the phase synchronization relations. Instead of single phase synchronization 
relations we have a phase synchronization space spanned by /3 which is the row space of n 
(33]. Testing for the rank of n can be viewed as testing for certain subspaces of the phase 
synchronization space represented by linear restrictions. 

The case r = d is a special and difficult one. In this case the process <f> t = {4>i > ■ • • > 4^)' 
no longer is cointegrated and there does not exist any equilibrium around which the process 
fluctuates in the above sense. Here some understanding is provided by the representation 
(10) which means that the d-dimensional process is driven by d-r integrated components. 
Thus in the case r = d there do not remain any stochastic trend terms who are driving the 
system and (loosely speaking) there is nothing left which can be synchronized. In fact the 
process in ® is stationary in this case. This is contradictory to our prior knowledge of the 
system: We know a priori that the phases constantly increase, i.e. they must have a trend. 
However the term ^ in ® does not produce any trend in case rank(n) = d. 

In that case there may perhaps be some other form of phase synchronization present in 
the data which we would call external synchronization caused by an external "pacemaker". 
Mathematically we mean by that two or more trend stationary processes, i.e. processes 



16 



with some intercept 7 and trend 5t plus some stationary part (regarded as noise in that 
situation). The phases 4>t anc ' <n tnen were externally m-.n- synchronized if m8\ = n5 2 . 
By using the model |[2) also oscillations different from pure cosine oscillations are possible. 
A trivial example for this are all daily oscillations which are caused by the daily cycle as an 
external pacemaker. Dependent on the specific situation, there may however exist examples 
of higher relevance. We emphasize that the detection of rank(n) = d is not sufficient for 
concluding to this form of external synchronization. Instead additional investigations are 
necessary which are beyond the topic of this article. For example one may fit a VAR-model 
with intercept and trend and test the relation m5i = n8i. 

To determine r = rank(n) we recommend using Johansen's LR-test described above 
based on the VEC-representation ©. It is important to include the parameter /j since this 
generates the deterministic trend oj and the phase shift (3 (as described above). If /3 is 
known one may use other testing methods such as the ADF-test. 

Kammerdiner and Pardalos (38l |37j also had the idea to use cointegration for phase 
synchronization and applied this to absence epilepsy data and to the analysis of neural 
data collected from primates during sensory-motor experiments. However their approach is 
different to the one presented here since they apply cointegration to the wrapped phases 
instead of the unwrapped phases as in our paper. We regard this as not adequate since 
the wrapped phases can hardly be regarded as realizations of an integrated process. As a 
consequence the phase processes in their examples were often tested to be stationary and 
the rank test for cointegration did often lead to the maximal rank r = d indicating stationarity 
of the phase processes instead of cointegration (of integrated processes). Probably as a 
consequence of these empirical findings the authors included the maximal rank r = d into 
their definition of phase synchronization and called this case perfectly synchronous. 

It is standard in physics to proceed as if the phases were observed directly and to ignore 
the effect of estimating the phases. The estimation of the phases is usually done by means 
of the Hilbert-transform on a segment. We will proceed in the following (and in particular 
in the example in Section [5) in the same way. In case of noisy systems a more realistic 
model may be to apply a nonlinear state space model as discussed in Appendix 1 where 
the unobserved phases form the state vector. However, the computations in that model are 
more challenging and time consuming. 
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Applying the method to phase processes 

We now summarize the main steps of the method. 

1 . Determination of the phases: 

Given the original observations of the oscillators, compute the phases processes <j>^ 
(J = 1, . . . ,d) by using the Hilbert transform or another method (see Section I2TT). Unwrap 
the phases prior to the subsequent analysis (that is the phases should not be restricted to 
[0, 27t)). For specific models such as 0-© other methods may be used. 

2. Testing for the rank of EE: 

Use the VAR-model as in (8) with /x ^ and conduct the Johansen LR-test as described 
above. We prefer the trace test in the "top -> bottom" - version. For the choice of the order 
p see the discussion in Appendix 2. It is also possible to test specific ft or hypothesis about 
p (35], Chapter 7). 

2a. Testing for cointegration when p is known: 

In the situation where the possible phase synchronization relations f3'4> t are known one may 
use instead the ADF-test (with a ^ and b = 0) to test for p'fc the null hypothesis integrated 
process vs. stationary process. If the test rejects we conclude to phase synchronization. 

This may in particular be applied for d = 2. In principle the use of our knowledge about 
P results in a slightly higher power of the test. However, our experience with the example in 
Section [5] indicates that the LR-test does not perform worse in this situation. 

3. Estimation, interpretation and further testing of the model: 

If not already provided from step 2, fit the model ® with (3 from step 2 to the data and esti- 
mate all coefficients. Check for significance of all coefficients and delete those coefficients 
which are not significant. Calculate uj and p as described above and write the whole model 
in the form 

p-i 

(AX t -oj ) = a (/3'X t _i -Po)+Yl T i ( AX *~i ~ "o) + Vt- 

with the estimated parameters - see ([18). Remember that the parameters a and ft are not 
uniquely determined (only the cointegration space is unique). This means that different /3 
can be chosen to formulate the same equilibrium relations (see also Juselius [36], Section 
8.5, "The cointegration rank: a difficult choice"). 

For testing the coefficients one can use the t-ratio. In the case of cointegration this t-ratio 
follows a standard t-distribution. We mention that the output of most statistical software 
automatically contains the t-ratio and the corresponding p-values. 
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3a. Testing for unidirectional coupling: 

An important special case of 3. is the testing for unidirectional coupling. If phase synchro- 
nization is detected one can test a for the direction of the dependency. For example in the 
VEC-model {14), where d = 2 the direction of coupling can be investigated with the 
adjustment coefficients o-i, a 2 . If, for example, a significance test suggests that a\ <0 and 
a 2 = one can conclude that influences (corrects) <$> but not the other way around. 
Thus, ffi is the driving force of phase synchronization and ffi is forced to adapt (unidirec- 
tional coupling). This means we have to test with the t-ratio whether a\, a 2 are significantly 
different from zero. 

For tests of more general hypothesis on the adjustment coefficients a see [35], Chapter 8; 
(361 . Chapter 11. 

Some software for performing these steps is discussed in the appendix. 

5 Example: A coupled Rossler-Lorenz system 

As an example we analyze with the above methods an unidirectionally coupled Rossler- 
Lorenz system. We proceed through steps 1. to 4. above. Computational details and 
additional modeling aspects can be found in Appendix 2. 

0. Simulation of the data: 

The system is defined through an autonomous Rossler attractor with configuration 

x\ = - 12(yi -Mi) + wi, 
y \ = 12(2:1 + 0.22/1), 

z\ = 12[0.2 + afi(a?i-5.7)] J 
and a driven Lorenz attractor given by 

x 2 = 16(y2 -x 2 ) — e(x 2 -xi) + w 2 , 
y 2 = 45.92x2 -2/2-^2^2, 
z 2 = x 2 y 2 -iz 2 . 

The parameters are taken from (27). The system noise variables w 1: w 2 are i.i.d. jV(0, 0.15 2 ) 
distributed similar to the system analyzed in [72J. The coupling is induced through the in- 
clusion of the xi term in the equation of x 2 with coupling strength e. Related systems have 
been discussed before (27] [68l [54l [56]. We have integrated the system with a step size of 
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Figure 1 : Rossler-Lorenz system: Trajectory of the Lorenz attractor in the xz-plane for coupling strengths 
e = (a) and e = 12 (b). (c) shows the trajectory in (b) corrupted by noise. 

0. 01 using a fourth-order Runge-Kutta algorithm [65]. After a burn-in period of 100 we have 
recorded 50 periods (5,000 data points) for the subsequent analysis. Different coupling 
strengths e are used (see below). 

To make the setting more realistic stationary observation noise is added to the oscillators 
before the phases are estimated: To the computed values Xi it ,yi,t,Zi,t at time t we add 
u\ = 0.9i4_i + v\ with v\ ~ Af(0, 0.1). 

The trajectory of the Lorenz attractor (in the xz-plane) for e = (no coupling) and e = 12 
(coupling) are given in Figure [Q It can be observed that, in contrast to the Rossler attractor 
(see Figure [2] (a)-(c)), the Lorenz attractor has two rotation centers in the uncoupled case. 
However, as a result of the symmetry the phase of the Lorenz attractor can be defined in the 
uz-p\ane where u = \Jx 2 + y 2 [64]. In the uz-plane only one rotation center exists which is 
illustrated in Figure (d)-(f). 

1 . Determination of the phases: 

We apply the Hilbert transform to estimate the phases as described in Section [2/Q The 
Hilbert transform is applied to the x and u coordinate of the noisy Rossler attractor and 
Lorenz attractor, respectively, using a rolling window of 1 ,000 data points. After the phase 
estimation we removed the first and last 500 data points (thus yielding 4,000 data points for 
each coordinate). 

Figure |3](a),(b) show 500 data points of the x and u coordinate of the Rossler and Lorenz 
attractors, respectively, for the uncoupled case (e = 0) and the coupled case (e = 12). The 
data are corrupted by noise as described above. In (c),(d) the corresponding (wrapped) 
phase estimates are given. It can be seen that in the uncoupled case the mean frequency 
of the Rossler system is smaller than the mean frequency of the Lorenz system. In fact, for 
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Figure 2: Rossler-Lorenz system: Trajectories of the Rossler attractor in the a^-plane (a)-(c) and Lorenz attrac- 
tor in the uz-plane (d)-(f) for coupling strengths e = (a),(d) and e = 12 (b),(e) (as a result of the unidirectional 
coupling (a) and (b) are identical). (c)/(f) show the trajectories in (b)/(e) corrupted by noise. 
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the given configuration the natural frequencies of the Rossler system and the Lorenz system 
are oj r = 0.129 and wl = 0.137, respectively (see below), i.e. the natural frequencies differ 
significantly in the uncoupled system. In (e),(f) the unwrapped phases are plotted and in 
(g),(h) the "stochastic trend" (in terms of cointegration), i.e. the unwrapped phases minus 
0.1292 are plotted (the latter being the deterministic trend of the Rossler-attractor and the 
coupled system as determined below). The shape of the stochastic trend in (g) already 
indicates no cointegration (phase synchronization) while in (h) it indicates cointegration. The 
oscillation of the phase of the Lorenz attractor is due to the fact that the Lorenz attractor has 
a slightly steeper ascent than a descent (in Appendix 2 we take this into account by choosing 
a higher model order of the AR-part). 

2. and 2a. Testing for phase synchronization: 

We first apply the Johansen rank test in the "top -> bottom" - version where we set the AR- 
order for simplicity to p = 2 (see Appendix 2 for a discussion). For example for e = 12 the 
hypothesis is clearly rejected with a LR-value of 92.44 (5%-critical value: 15.41) and the 
hypothesis H(l) is clearly not rejected (i.e. accepted) with a LR-value of 0.01 (5%-critical 
value: 3.84). Thus the test reveals rankll = 1, and we conclude to phase synchronization. 
Furthermore, the procedure also detects the 1 : 1 - relation of the phase synchronization from 
the corresponding eigenvector. 

FigureH(a) shows the values of the H(0) - test statistic for different e. The dashed line and 
the dotted line are the critical values 15.41 and 19.62 for the significance levels 0.05 and 0.01, 
respectively [36], Appendix A, Case 3). For those e where the test rejects the hypothesis 
H(l) was tested afterwards and in all cases not rejected. Thus the plotted value of the T-L(0) 
- statistic determined the decision phase synchronization vs. no phase synchronization. The 
values for e > 11.1 are all highly significant corresponding to phase synchronization while 
the values for e < 11.1 are nearly all not significant indicating no phase synchronization. 

In Figure H](c) the phase synchronization index R 2 as a classical method is plotted for 
the same e. The associated critical values are computed as explained in [72j (this turned 
out to be difficult - see Appendix 2). Basically the results of the two tests agree in all cases 
apart from e = 9.6, 10.2, 11.0, 11.1 (and 10.6 as a "borderline"-case). The difference can in 
all 4 cases be explained by 2n -jumps of the phase - differences which are penalized by the 
LR-test but not penalized by the i? 2 -statistic. The situation is highlighted in Figure [5] where 
we have plotted the estimated phase differences (fry—^y for e equal to 0, 9.6, 10.2 and 12 
(the situation for 11.0, 11.1 being similar to the cases 9.6 and 10.2). We think it is a matter of 
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Figure 3: Rossler-Lorenz system: 500 data points of the x coordinate of the Rossler attractor (solid lines) and 
of the u coordinate of Lorenz attractor (dashed lines) for e = (a) and e = 12 (b). The data points are corrupted 
by noise, (c), (d) give the corresponding wrapped phase estimates, (e), (f) the unwrapped phases, and (g), (h) 
the stochastic trend of the unwrapped phases. 
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Figure 4: Rossler-Lorenz system: (a) values of the LR-test for H(0), (b) values of the ADF-test, and (c) the 
phase synchronization index R 2 , for coupling strengths e between 8 and 13. The dashed lines and the dotted 
lines correspond to the significance levels a = 0.05 and a = 0.01, respectively. 
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Figure 5: Rossler-Lorenz system: Phase differences 4>t ~ 4>t for coupling strength e = (black), e = 9.6 
(green), e = 10.2 (red) and e = 12 (blue). 

taste whether one calls the cases 9.6, 10.2 phase synchronized. 

We have also used the testing procedure as described under 2a. based on the ADF-test. 
Here we have to specify the possible phase synchronization relation in advance. Based on 
the plot in Figure [3] (b) we have decided to test for a 1 : 1 synchronization relation, that is 
we have applied the ADF-test to </>f'—</>i 2 ' as described above. The results for different e 
are given in Figure H(b) with the vertical axis being upside down. The dashed line and the 
dotted line are the critical values -1.65 and -2.33 associated with the significance levels 
0.05 and 0.01, respectively (from the standard normal distribution - see |28], p. 529, Case 3). 
The results strongly coincide with the LR-test (with a slight difference for e = 10.5, 10.6, 10.7). 

It is remarkable that the LR-test and the ADF-test as two standard tests from the theory 
of cointegration perform as well as the i? 2 -test which has been explicitly tailored for phase 
synchronization. The outcome is only different for 4 out of 51 values and this behavior can 
be explained in all cases by the 27r-jump penalization. We mention that this difference may 
even be cured by including structural shifts of size 2ir into the cointegration approach similar 
to |48]. 

Since the LR-test seems to be a little bit closer to the i? 2 -test and its setting is more gen- 
eral (no pre-specification of the synchronization relation; arbitrary dimension, i.e. networks) 
we have the tendency to prefer the LR-test. 

3. and 3a. Estimation, further testing, and unidirectional coupling: 

Given that 4>f \ are cointegrated with cointegrating vector j3 = (1, -1)' we estimate the 
VEC-model (13}, ([15) to analyze the joint dynamics with the aim to uncover the directional 
coupling. Here we restrict to the case e = 12. The estimation results are given in TableQ] 
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Parameter 


Estimate 


Std. Error 


t- ratio 


p-value 


cy-t 


0009 


0.0015 


0.62 


0.53 


n 


0.059 


0.016 


3.71 


0002 


8-1 


0.024 


0.014 


1.70 


0.089 


1 1 1 


0.1177 


0029 


41.12 


< 0001 




0.0110 


0.0011 


9.65 


< 0.0001 


72 


0.726 


0.011 


68.38 


< 0.0001 




0.028 


0.012 


2.38 


0.017 


M2 


0.0262 


0.0022 


12.11 


< 0.0001 



Table 1 : Estimated VEC-model for the Rossler-Lorenz system. 

In the equation of A^ x) only m and 71 are having a significant t-ratio. In the equation 
of &<j>y all parameters except 5 2 are highly significant. This result agrees strongly with 
the properties of the investigated Rossler-Lorenz system. First, the VEC-model reliably 
identifies the direction of the coupling because «i = and a 2 >0. Second, 5i and 5 2 as well 
as the cross-correlations of the residuals with Corr (77^, rj^ ) = 0.022 are not significantly 
different from zero. This indicates that the dependency of the phases is purely explained 
by the error correction mechanism. This agrees with the type of coupling of the oscillators. 
Finally, a reestimation of the model with only the significant parameters yields 



(i) 



0.060 A^i\ +0.121 + 



(2) 



0.011 



+ 0.727 + 0.030 + i)f ] . 



(16) 
(17) 



Note also the more informative representation {T8]) below and its discussion. 

We now discuss the estimated system in terms of the modified VEC-representation (13) 
and the MA-representation |T7j|. At the same time this system is a nice illustration of the 
theoretical setting presented in Section [3J] (and a nice example for a cointegrated system). 
With the above estimates we obtain 



1 



Pi 



a 





0.011 



0.121 
0.030 



This leads with straightforward calculation to the other values 



Ti 



0.06 
0.73 



0.94 
0.27 



C = 0.94" 



'1,0), uj = Ch = 0.129 



and A) = 0.49. With these values we obtain the VEC-representation 



A<A, 

A<A 



0.129 
0.129 





0.011 



-0.49 + 



0.06 
0.73 



A#\ 



0.129 
0.129 



(18) 
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and the MA-representation 

= (1) X>* W + (aisg) 1 + C ^ L ^ + /*) + (19) 

(ED has the meaning that the phase process is pulled towards the equilibrium defined by 

P'4>t-PQ = 4>t -<t>t - °- 49 = w ' tn tne f° rce a = (oon) wn ' cn 9 ets active as soon as the 
process leaves the equilibrium. In this case a\ = i.e. the correction is on process 2 only. 
This means that the unidirectional coupling between the Rossler and the Lorenz system has 
been detected correctly. The detection was implicitly based on a test since the coefficient a\ 
in TableOQwas not significantly different from 0. 

The equilibrium </>| 1) -<^ 2) - 0.49 = means that the phase of the Rossler system is ahead 
of the phase of the Lorenz system in the average by 0.49. This is nicely confirmed by the 
plot in Figure |3](d). 

The representation 09) has the meaning that both phases are mainly pushed in direction 
(I) - both with a deterministic and a random walk component. The stationary component 
C*(L)(r] t +fi) (whose precise form is complicated) acts as a disturbance component pushing 
the system with little shocks again and again out of the equilibrium. 

For comparison, we also fit the phase model A<p t = -yAfo-i + n + rjt to each oscillator in 
the uncoupled case e = which gives 

(A^ 1} - 0.129) = 0.060 (A^i\ - 0.129) + $\ (20) 
(A4 2) - 0.137) = 0.735 (A^\ - 0.137) + r,f\ (21) 

( ([16) and <[20j> are of course identical). As mentioned above the natural frequencies of the 
Rossler and the Lorenz system therefore are ojr = 0.129 and vol = 0.137, respectively, 
i.e. the natural frequencies differ significantly in the uncoupled system, while in the coupled 
system they both are u = 0.129 (remember that ut is the deterministic part of the trend). 

One effect should be kept in mind: The prior application of the Hilbert transform leads to 
phases which are smoother than the original "true" ones. This effect needs to be investigated 
in the future (see Appendix 1). As a consequence the estimated variances of r) t and & 
will usually be smaller than the true ones. In <(20]>, ((21} one finds Var^ 1 ^ = 0.0034 and 
Varr/ t (2) = 0.0025 leading to Var (A^ 1} )= 0.0034 and Var (A^ 2) )= 0.0054. This means that 
the Lorenz attractor has a more varying frequency than the Rossler attractor and it is likely 
that this also holds for the true frequencies. 
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6 Conclusion and outlook 



We have pointed out the connection between the theory of cointegration and the theory 
of phase synchronization. In particular a cointegrated dynamical system can be used as 
a stochastic model for a multivariate phase process which describes the behavior of the 
phases in detail. Contrary to other concepts like the spectral coherence which only describe 
the properties of the phase processes, we now have a specific model in the time domain for 
the dynamics of the phase processes. The model leads for example to the identification of 
the equilibrium relations (phase synchronization relations) from data including phase shifts, 
to the identification of the coupling structure (by using tests for unidirectional and bidirec- 
tional coupling), to the identification of phase synchronization spaces in networks, and to 
tests for phase synchronization itself. We have demonstrated this for a coupled Rossler- 
Lorenz system with noise. 

In particular in neuroscience where phase synchronization is regarded as essential for 
functional coupling of different brain areas the new methods coming from cointegration may 
enhance the present methods like spectral analysis, correlation and coherence analysis, 
triplet analysis, joint phase histograms etc. This however needs to be investigated and 
confirmed. 

From a physical view the potential of cointegration for phase synchronization can perhaps 
best be understood from a comparison with the Kuramoto model for phase synchronization. 
As pointed out in this paper cointegration can be understood as a stochastic Kuramoto-type 
model which can be fitted to the data. The cointegrated system describes both: 1) how the 
system is constantly pulled towards the equilibrium determined by the phase synchronization 
relations; and 2) how the system is constantly disturbed by little shocks and kicked out out 
of the equilibrium. It seems to be interesting to reinvestigate the properties of the classical 
Kuramoto model in this new framework. 

Another important point is that the cointegrated phase model covers both: the case where 
the phases are synchronized and the case where the phases are not synchronized. As a 
consequence the diagnoses on this can be made from a statistical fit of the model to the 
data. In networks cointegration allows for different types of phase coupling which can also 
be diagnosed from data. 

At the same time the comparison to the Kuramoto model shows the difference to previous 
research: In the classical view (say when the phase synchronization index R 2 is calculated) 
phase differences of 2tt are not penalized, that is two processes whose phase difference 
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moves (quickly) by 2ir may be regarded as phase synchronized. On the contrary such 
processes would not be recognized as phase synchronized in the cointegration approach of 
this paper. This would require the inclusion of structural shifts into the cointegration model 
which is possible. In our simulations this difference occurred in 4 out of 51 cases and a 
closer inspection revealed that it is not clear at all what the better method is in this situation. 
Moreover, we guess that this situation hardly occurs in most applied data sets. 

In the theory of phase synchronization many topics have been investigated which have not 
been looked at in cointegration. Some of those topics may be of interest in economics and 
worth to be looked at in the different setting of cointegration. For example in the framework 
of the Kuramoto model the detailed investigations on the conditions under which oscillators 
are able to synchronize, on the onset of synchronization, on the mean field, or the thermo- 
dynamic limit when the number of oscillators tends to infinity ([3], Chapter III). In that context 
it may for example be of interest to study cointegrated systems where the dimension d tends 
to infinity. 

7 Appendix 1 : VEC - state oscillators 

Inspired by the motivation at the end of the introduction we define the class of 
VEC - state oscillators as a general model for (possibly phase synchronized) d-dimensional 
oscillators with random phases. The model is the nonlinear state space model with 

State Equation: 

A<p t = n <p t -i + r * A( ^- +P + Vt> Vt- AA(0, E„) (22) 
i=i 

Observation Equation: 

Y t {l) = a { co S (4 l) ) +bi + ef, e t ~ Af(0, S e ) (23) 

where e t and rj t are mutually and serially independent. In its simplest form the amplitude a 
and the baseline b are parameter vectors. In some cases they need to be time varying (see 
(ii) below). The case rank(n) = is included (the case rank(n) = d requires the additional 
term vt in ((22) for a meaningful phase-model - see Section 13}. If 1 < rank(n) < d - 1 we 
have phase synchronization. In that case we prefer writing the state equation in the form 

p-i 

(A& - oj ) = a (P'<p t -i ~Po)+Yl T i i^t-i - u o) + Vt (24) 
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(see Section |3) with the phase synchronization relations p'fa —fa = 0, the mean phase shift 
fa, and the deterministic trend uot. n is now decomposed into n = a/3' with dxr - matrices 
a and /3 of full rank r. 

The most common approach in physics is to estimate the phase process 4>f from the 
oscillator by means of the Hilbert transform. In terms of state space models one may 
regard this as a smoother for estimating the mean of the conditional distribution (although 
the properties are not clear). The approach of this paper is to use the estimated phases "as 
if they were observed" and to apply standard cointegration techniques to it. 

Of course the effect of the Hilbert transform on the properties of these estimates is not 
clear, and the significance levels of the tests need to be investigated. This is beyond the 
scope of this paper, but we want to provide some heuristic arguments why everything re- 
mains the same even in this situation: Taking the conditional expectation (given Yi, . . . ,Y n ) 
in ([24) reveals that the linear structure of the conditional means is the same as before. The 
difference is that the variance of the innovations will be smaller (the estimated phases are 
smoother than the original ones) and the errors become dependent. A least squares re- 
gression based on the conditional expectations should also lead to consistent estimates for 
the parameters and in particular the t-statistic for a specific parameter (which is rescaled 
by the empirical variance) should follow asymptotically a t-distribution. The argument is the 
same for the ADF-test in the case a / 0, b = since this is the "regular" case where the 
drift dominates the stochastic trend and the asymptotic distribution of the unit root estimate 
is Gaussian ([23], Theorem 10.1 .5). 

For the LR-test with its nonstandard distribution a different heuristics is needed: Inspection 
of the proof of the asymptotic distribution of the LR-test ([35], Theorem 11.1) reveals that its 
limiting distribution only depends on the number d—r of common stochastic trends and on 
the model for the deterministic terms. The model for the deterministic terms is fixed - so we 
only have to check that the number of stochastic trends remains the same for the estimated 
phases (estimated with the Hilbert transform) as for the original unobserved phases. To 
check this we consider in this heuristics only the case of 1 : 1 phase synchronization relations, 
i.e. where the system splits up into d-r oscillators and each of the remaining d oscillators 
is stuck to one of the former ones by the cointegration relation, i.e. we have d-r groups 
of oscillators of different size (this is a common case for phase synchronization). We now 
have only to check that the groups of phase processes "stay together" under the Hilbert 
transform. First it is obvious that an integrated process stays integrated: The phases will 
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constantly increase and a trend stationary process can be excluded because this is the 
case where the length of the cycle is almost stuck to a pregiven length and it is clear that a 
true phase process with variable cycles will not be transformed by the Hilbert transform into 
a phase process with almost constant cycle-length. In addition it is also clear that the phase 
processes will stay in the same group, i.e. they will not jump to another group with a different 
trend since this is clearly reflected by the number of cycles. Thus the number of stochastic 
trend should stay the same under the Hilbert transform. 

We emphasize that these are only heuristic arguments. A mathematical proof that the 
asymptotic distributions of the test statistics stay the same would be highly welcome. It is 
obvious that deriving such a result may be very challenging. 

A more sophisticated approach seems to be to estimate the phases in the above model 
by means of a particle filter. This has been done in [12] in the univariate case d = 1 (in the 
more general setting (i)-(iii) from below). However, also with this approach the situation is 
not clear: A proper test for synchronization would then e.g. be a likelihood ratio test based 
on the original observations Y t . The test statistic for such a test could be approximated by 
means of a particle filter - but its distribution under the null hypothesis is also not clear and 
difficult to derive. 

Our personal view is that for systems with stronger noise the above particle filter or a 
periodogram based method lead to better estimation results. Here it is worthwhile to study 
the questions raised above. 

For chaotic oscillators as in Section [5]the Hilbert transform seems to be the better choice 
(furthermore it is used in most applications). Here the misspecification error has to be taken 
into account when the above questions are investigated (see Appendix 2 below). 

Generalizations of the state space model 

(i) Other oscillation patterns: 

For non-cosine type signals one may use instead the observation equation 

Y t = a t f(4>t) + h + e t , 

with / being a 27r-periodic real valued function representing the oscillation pattern [T2]. / 
typically is unknown and needs to be estimated from the data. An example are ECG-data. 

(ii) Modeling the baseline and the amplitude: 

Dependent on the data one may need a time varying amplitude and a baseline. This can be 
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achieved with a state space model also for a t and b t , for example a VAR-model with mean 
different from zero. 

(iii) Positivity of the phase increments: 

In nearly all cases it is clear that the phase increments of an oscillator should be positive. 
However, in the above model with Gaussian noise also negative phase increments may 
occur. This happens with small probability since w usually is large in comparison with the 
standard deviation of the r) t . Furthermore, this should not be a big problem if the model is 
mainly used for estimation and testing (and not for simulation). 

To overcome this problem we have used in [12] in the case of a single oscillator an in- 
tegrated ACD (autoregressive conditional duration) model for the state. At present we are 
having no idea how this ACD-model can be extended to the cointegrated case (see however 
[52] for bivariate ACD-models). 

In [T2] all three generalizations have been combined in the univariate case. The phases 
have in that case been estimated by a fixed-lag particle smoother and the unknown oscilla- 
tion pattern / by a nonparametric EM-algorithm. 

8 Appendix 2: Computational and modeling aspects 
8.1 Computational aspects 

There exists several software programmes on cointegration analysis. For example the CATS 
module within the RATS Econometrics Software ( [http://www.estima.com/} and OxMetrics 
(http://www.oxmetrics.net/). For the analysis of the Rossler-Lorenz system in Section [5]we 
have used the R-packages urea (http://cran.r-project.org/web/packages/urca/index.html - 
see [60]) and systemfit (http://cran.r-project.org/web/packages/systemfit/index.html - see 
[31]). Here are some details: 

1. Initial situation: The unwrapped phases are given as vectors phil, phi2 (e.g. determined 
by the Hilbert-transform). The first values phi 1 [1 ] and phi2[2] should both lie in [0, 2w). 

2. Johansen LR-test (top-bottom version): 
require (urea) 

joha.test <- ca.jo(cbind(phi1 ,phi2),type = "trace", ecdet = "none", K = 2, spec="transitory") 
print(summary(joha.test)) 

Note: The specification ecdet = "none" fits © with ^ i.e. with w and /3 (Case 3 in [36] . 
Chapter 6.3) (while ecdet = "const" and ecdet = "trend" correspond to case 2 and case 4, 
respectively). We recommend using the quantiles from the appendix of [36] (Case 3) instead 
of the printed values in the summary. 
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2a. ADF-test for testing phase synchronization 

(testing the phase-difference for integrated process vs. stationary process): 
require (urea) 

beta2 <- 1 # Specify cointegration vector beta 
adf.test <- ur.df(phi1 - beta2 * phi2, type="drift") 
print(summary(adf.test)) 

Note: The specification type="drift" corresponds to Case 3 in (28), p. 529. The relevant 
statistics is tau2 which should be compared with the quantiles of the standard normal distri- 
bution (for sample sizes less than 300 the t-distribution probably gives the better quantiles). 
The quantiles printed in the summary correspond to Case 2 in [28] where a different hypoth- 
esis is tested with the same statistics. 

3. and 3a. Estimation, further testing, and unidirectional coupling: 
require(systemfit) 

T <- length(phM) 

err.lagl <- (phi1-beta2*phi2)[-c(T-1 , T)] 
err.lag2 <- (phi1-beta2*phi2)[-c(T-1 , T)] 
dphil <-diff(phi1) 
dphi2 <- diff(phi2) 

diff.dat <- data.frame(embed(cbind(dphi1 , dphi2), 2)) 
colnames(diff.dat) <- c('dphi1', 'dphi2', 'dphil. 1', 'dphi2.1') 
eqPhil <- dphil ~ err.lagl + dphil. 1 +dphi2.1 
eqPhi2 <- dphi2 ~ err.lag2 + dphi2.1 + dphil .1 
system <- Iist(phi1 = eqPhil, phi2 = eqPhi2) 
estSystem <- systemfit(system, data=diff.dat) 
print(summary(estSystem)) 

4. Final estimation of the reduced system: 
eqPhil <- dphil ~ dphil. 1 

eqPhi2 <- dphi2 ~ err.lag2 + dphi2.1 
system <- Iist(phi1 = eqPhil, phi2 = eqPhi2) 
estReducedSystem <- systemfit(system, data=diff.dat) 
print(summary(estReducedSystem)) 

Finally, we mention that the critical values for R 2 are not tabulated. They depend on certain 
parameters which need to be estimated based on the phases [72]. In our example, we found 
that the critical values are sensitive to the segment length parameter and we struggled to 
obtain the critical values plotted in Figure 0(c); we still remain a bit unsure about the plotted 
values but they seem to be reasonable. 

8.2 Modeling aspects 

One needs to be aware of the fact that the VEC-model has been applied in Section [5]in the 
highly misspecified situation of two noisy chaotic oscillators. Our goal was to demonstrate 
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that already a "simple" VEC-model can be used to successfully identify phase synchroniza- 
tion and unidirectional coupling. For that reason we have chosen the fixed order p = 2. 

In a more sophisticated modeling one needs to select the order p. A common suggestion 
is to use a model selection criterion for choosing p such as the AIC or the BIC (cf. (25] , 
p.752). For phase synchronization we clearly prefer the BIC instead of the AIC since the 
AIC since is not consistent and the model order is typically too large for the present purpose. 
Without having tested this, we are convinced that the BIC works reasonably well in many 
situations of phase synchronization. 

However, the situation is different for the Rossler-Lorenz system investigated in Section 
To understand the problem one has to realize that all model selection criteria have been 
developed under the premise that the true model lies within or at least close to the fitted 
model class. In the present highly misspecified situation this is not true. Accordingly we 
have found with the BIC the exorbitant model order p = 31 and with the AIC even p = 95 for 
the phases of the Lorenz attractor (although not all lags were significant). The result means 
that we need a very high order to obtain a good approximation of the misspecified system 
with a VEC-model. On the other hand we would hardly trust the outcome of significance 
tests etc. with such a larger order. For that reason one has to refrain from these automatic 
procedures for chaotic oscillators (or has to develop new ones which perform better in such 
situations). 

There is a remarkable feature in Figure [2] (g),(h), namely the cosine-type oscillation in 
the stochastic trend of the phase of the Lorenz attractor. This oscillation is due to the fact 
that the Lorenz attractor has a slightly steeper ascent than a descent. From a time series 
perspective it is clear that this can be modeled by a higher order p of the AR-part - the rule of 
thumb saying that each cosine cycle can be modeled properly by 2 extra model orders of p. 
This of course is ad hoc and we have fitted AR(p)-models with different model orders where 
indeed the increase from p = 2 to p = 4 showed a remarkable improvement of the fit for the 
Lorenz attractor while higher model orders than p = 4 only showed smaller improvements. 

We therefore have redone the complete analysis of Section [5]with this order. There was 
a clear quantitative improvement but qualitatively all results stayed the same: In particular 

- the outcome of all tests for phase synchronization was essentially the same (although 

some decisions could be made with higher significance); qualitatively the cases e = 
9.6, 10.2, 11.0, 11.1 discussed above remained different to the i? 2 -statistic; 

- the findings on the undirectional coupling remained the same - that is ai was again tested 

to be while a 2 was clearly different from leading to unidirectional coupling; 
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- we obtained again w = (o 129) anc ' ^0 = 0.48, i.e. the same deterministic trend and the 
same phase synchronization relation. 



The model equation was 



-w 





0.018 



+ 




0.21 



(i) 

t-2 



t-2 



0.48 + 
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t-i 
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-CJ 



-w 



Vt 



meaning that all structural aspects of this model are the same as for p = 2. 

Another consequence of the highly misspecified situation is that it makes no sense to test 
at the beginning the hypothesis integrated process vs. trend stationarity with the ADF-test 
or the hypothesis trend stationarity vs. integrated process with the KPSS-test. We have 
tried this and found many situations where both tests reject. Usually this is contradictory. 
In the present situation it just means that the true process is neither trend stationary nor an 
integrated process but the phases of an chaotic oscillator with some noise. 

The consequence is that one has to make the decision whether it makes sense to apply 
cointergration from common sense, meaning that one has to exclude those cases where the 
length of the cycle is almost stuck to a pregiven length (for example all daily cycles have to 
be excluded since one knows long in advance where the phase will be at a certain time). 
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